
%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------


mhrund    = 1; 


%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;

% parameter names
pnames = feval(fnames_.pnames,msel);
mhrun = 1;
mhrunid  = num2str(mhrun','%2.0f');

block1   = 11;		% starting block
blockn   = 150;		% ending block
nsim     = 10000;	% number of simulation draws in each block

%--------------------------------------------------------------------
% means and standard deviations
%--------------------------------------------------------------------


 lpath = [path_.para 'mhrunpm' mhrunid  '/' ];

	% load file
	lfile = [path_.para 'stat_mhrunch01.mat'];
	load(lfile);

%alpha1_grid = 0.85:0.05:2.5; 
%sd_f_grid   = 0.0015:0.0001:0.0035; 
p_grid = 0.5:0.01:0.99;
lk = length(p_grid); 
para = drawmean';
wu1 = zeros(lk,100*(blockn-block1+1)); 
wu2 = zeros(lk,100*(blockn-block1+1)); 
wuint = zeros(2,lk); 
wuint2 = zeros(2,lk); 
k=1; 
  for jblock=block1:blockn

	% load file
	lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
	load(lfile);


    % number of simulations in each block
	
   for jj=1:100:nsim
        para=parasim(jj,:);  
for j=1:lk 
%alpha1 = alpha1_grid(j);
alpha1 = para(1); 
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);    
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a  = para(8);
rho_f  = para(9);
rho_i  = para(10);
sd_a1   = para(11);
sd_f1   = para(12);
sd_i1   = para(13);
sd_a2   = para(14);
sd_f2   = para(15);
sd_i2   = para(16);
piess  = para(17);
gy     = para(18);
lnA_0  = para(19);
yst    = para(20);
p11 = p_grid(j);
p22 = para(22);
p11_f = para(23); 
p22_f = para(24); 



P_i = [p11 1-p11;
       1-p22 p22];
   
P_f = [p11_f 1-p11_f;
       1-p22_f p22_f];   

P = kron(P_i,P_f);

iss    = gy-log(beta1) + piess;

%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------

  A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
        -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];

    solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

    A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
        -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];
    solvec_f = inv(A_Af_mat) * [0;0;-kappa1;-kappa1;0;0];

    A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
        -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];
    solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];

   
    pi1a = solvec_a(3);
    pi1f = solvec_f(3);
    pi1i = solvec_i(3);
   

    var_a1 = (sd_a1^2)/(1-rho_a^2);
    var_u1 = (sd_f1^2)/(1-rho_f^2);
    var_e1 = (sd_i1^2)/(1-rho_i^2);


    wu1(j,k)  =((pi1f^2)*var_u1)/(((pi1a^2)*var_a1)+((pi1f^2)*var_u1)+((pi1i^2)*var_e1));
    wu2(j,k)  =((pi1i^2)*var_e1)/(((pi1a^2)*var_a1)+((pi1f^2)*var_u1)+((pi1i^2)*var_e1));
end
    k=k+1; 
   end
  end 

for jk=1:lk 
   wuint(:,jk)=hpdint(wu1(jk,:)',0.9);
   wuint2(:,jk) = hpdint(wu2(jk,:)',0.9); 
end 
  linewidth=2;
linecolor1='r';
linestyle1='-';
linestyle2='--';
xxtp1=   1.79599541;  
xxtp2=   1.14322405; 
ymin1=0; ymax1=1;
ymin2=0; ymax2=1; 


subplot(2,1,1); plot(p_grid,[mean(wu1,2) wuint(1,:)' wuint(2,:)'],'k',[0.5 0.5], [0 1],'--k',[1 1],[0 1],'--k'); 
%line([xxtp2 xxtp2],[ymin2 ymax2],'LineStyle',linestyle2,'Color',linecolor1,'LineWidth',1);
line([xxtp1 xxtp1],[ymin1 ymax1],'LineStyle',linestyle1,'Color',linecolor1,'LineWidth',2);
subplot(2,1,2); plot(p_grid,[mean(wu2,2) wuint2(1,:)' wuint2(2,:)'],'k',[0.5 0.5],[0 1],'--k',[1 1],[0 1],'--k'); 
%line([xxtp2 xxtp2],[ymin2 ymax2],'LineStyle',linestyle2,'Color',linecolor1,'LineWidth',1);
line([xxtp1 xxtp1],[ymin1 ymax1],'LineStyle',linestyle1,'Color',linecolor1,'LineWidth',2);


			